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ABSTRACT 


The analysis of nonlinear reactor kinetics with Doppler feedback 
4 
is developed by the use of the Poincare-Lindstedt perturbation technique. 
The two types of temperature dependent reactivity coefficients to be 


investigated are 


dk ] dk 1 
ar * q Be dT r. 


Each of these reactivity models are analyzed for the reactor with positive 
and negative coefficients of reactivity and infinite slab and finite cylin- 
drical geometries. The Poincare-Lindstedt perturbation solution is con- 
structed in the form of an Ben about the unperturbed state of the 
reactor. Each of the perturbed equilibrium solutions is subsequently 


analyzed for stability. 
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1. THE NONENNEAR TEMPERATURE DEPENDENT REACTOR 


jae GOVERNING EQUATIONS 
In a homogeneous, one-velocity, bare reactor the neutron popu- 


lation is described by the equation (1) 


1 oO 2 
E URI DP ext) EPA E 


N 
R (1-3) XP (x, t) + p ) I ES 


ne first term Eq. (1) represents the time rate of change of the neutron 
flux within a differential control of volume. This term is balanced by the 
net effect of the diffusion of neutrons, the absorption of neutrons, and 
the production of both prompt and delayed neutrons within the volume. 


The concentration of the precursors is described by [1] 


ORG, +) | 
J a E + R Laß: PX, 1) (2) 
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The solution of these equations with temperature dependent parameters 
represents the goal of this work. 
The temporal and spacial coordinates may be nondimensionalized 


by introducing the definitions 


= yn 


K D X A = Wee t ; 





Equations (1) and (2) become the following 
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where the precursor concentration and the flux must have a zero value 
on the dates of the reactor. The solution is desired for the infinite 
slab reactor and the finite cylindrical reactor. The first of these cases 
is the slab which may be envisioned as a reactor which is infinite in 
two dimensions and bounded by the coordinates +1 and -1 in the third 
dimension. This simple geometry will easily lend itself to analysis, 
but permits the observance of mathematical rigor and the physical inter- 
pretation of the results. The boundary conditions require that the flux 
is zero at the coordinates + 1. The second caSe is that described by 

a finite cylinder which has boundaries on the longitudinal axis at +l 
and at a nondimensional radial distance of To: Therefore, the flux is 
required to be zero at the extremes of the longitudinal coordinate, 


á= + l, and at the radius r = ro 


Br DOPPLER IEE Pease 
The lumped parameter k is generally referred to as the neutron 
multiplication factor. In elementary reactor analysis the multiplication 


factor is evaluated as a constant given by the product [1] 





2 = R Pr Y 

However, in reactors that operate at high temperatures the multiplication 
factor becomes temperature dependent as will subsequently be shown. 

One of the major considerations to be made when accounting for 
the temperature dependence of the multiplication factor is that described 
by the nuclear Doppler effect. Physically, this phenomenon results from 
the temperature broadening of the capture cross-sections of the fuel 
nucleii. Thus the Doppler effect may be envisioned by a temperature 
dependent competition for neutrons by capture and fission processes. 
The quantitative description of the Doppler effect is achieved by making 
the following assumptions: 

a. The fuel density as a function of velocity is given by a 
Maxwellian distribution function. 

Dr The capture cross-sections near an isolated resonance is 
given by the Breit-Wigner formula. 
It may be shown that the average radiative cross-section, O y : 
and scattering cross-section, Os , are described as a function of 


energy and temperature given by [1] 


oy Tr a 
Oy (E,T) = we W (t,x! 


and 


p A a 
A TE 444273 ‚a NR) + ATR 





where WY ( 81%) and Ie ( I) are tabulated functions of 





temperature, 2 
x = 7 (E-&) 
Po 
E AE, RT 
A 
The effectof Oy and O, is most conveniently incorporated 


into one quantity by introducing the effective resonance integral, I, 
defined by [2] 


o Je OO dE 
Il = Tie cee E One ©, 107. E 


where OÖ, is the additional scattering cross-section per fuel atom. 
It is usually more convenient to analyze the temperature dependence 
through the J function which is related to the resonance integral as 


follows [1] 





= Gr F(t) 
ae f 
where 5 Na Gon y 
2 Dre 
and the subscript M refers to the moderator. 
If attention is returned to the concept of the multiplication factor, 
its dependence upon temperature may be investigated [1] 


Since 
R « P 
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and 
+p = exp Ge) 
the Doppler effect on k is studied by considering the temperature effect 
on J . However, the preceding analysis showed that 
om Ty 


I = Er Jr (834) 


and it may be simply concluded that 


Ro exp [IB 
If the derivative with respect to temperature is taken, it is obvious that 


\ dk AI 


— Oo a (5) 
ee oon dT 


However, under most operating conditions the value of p is near 


unity. Thus for small changes in k the following relation is valid, 


dk „ As 
a T oul 


(Sa) 


The left hand side of Eq. (5) is referred to as the temperature coefficient 
of reactivity and provides a functional dependence upon the of function 
and, in turn, upon the changes in temperature of the reactor. 

In a paper by A. Reichel the temperature dependence of the 
function was investigated [2]. The results were based upon evidence 


2 
gathered from experimentation with U > and from analytic investigations 


21 





into the nature of the J function. Reichel concluded that physically 
realizable states of the reactor must be located in the range 


MOLES < © and MODO! <3< 10. 


The results are most simply presented by Figure 1. The variable 


5 
ß is described by K where K = Am UOR) £ 
km (2) 


The shaded region represents states which are physically realizable for 


yess. 


From the analytic investigations it was observed that values of 
ß at the extremes of its range coincided with conditions which 

required that the y function be a constant. Therefore, the tempera- 
ture dependence of Jd is destroyed and, in turn, the Doppler effect 
and temperature dependence of the coefficient of reactivity is nonexistant. 
In other regions of Figure 1 the dependence of % was found to vary 
with temperature raised to a constant power. These regions are indicated 
on Figure 1 by the lines labeled with Bu Thus for any particular state 
of operation defined by values of R and E the dependence of 
upon temperature may be defined. If Eq. (6) is evoked, this Tohe- 
pendence is transferred to the governing equations, Eqs. (3) and (4), 
via the temperature coefficient of reactivity. 

As may be noted from Figure 1 the variance of the J function 
with temperature may take on many different forms. However, in an 


analysis presented by Hummel and Okrent [3] the ig function was 


found to vary between the functions 


I «MT Toa y 
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These results are used to determine the functional dependence of R 


with temperature. The results are that 
-į 
IR e T 
(6) 


and 


a 


T Pr 


a 


1 
respectively. In the remainder of this paper the effect of these types 
of temperature feedback are analyzed. 
The proportionalities in Eq. (6) are removed by the addition of a 


constant. The values of are determined by experiments per- 


formed on the reactor in question. Some representative values are given 


for various reactors. 


TABLE 1 [4] 
Nuetron 
Reactor Lifetime Temperature Coefficient 
Fermi 14 A sec -3.4x No N k/k per °C 
EBR II .68 p sec Deo A k/k per °C 
Hallam .37 m sec + .3x 1079 A k/k per °F 
Peach Bottom .23 m sec -2.6 x or Aer © 
EGCR .40 m sec -4.5 X 10 NEP e 
EBOR e227 je sec -1.5 x re i pee © 


1 

In the subsequent analysis the need to differentiate betwecn the 
feedback models will arise. The distinction will be made by referring 
to the "T71 model" or mp =3/2 model." 


14 





3/2 


= = 
Bor the TI and the T model reactors the temperature dependence of 


the multiplication factor may be described by the equalities 


AR au 


aT 7 T 
and 4 o = (7) 
a 
ar “7 TW 


respectively. It may be noted that the negative sign is incorporated 

into above equation. This convention is generally observed because of 

the greater frequency of reactors with negative temperature coefficients. 
| ' > | AR 

In Eq. (7) the units of A and Q. are defined such that ar 


fers the units of A k/k per °C. 


©. ASSUMPTIONS 

Before the analysis of the reactor may be undertaken two further 
assumptions must be made concerning the dependence of the neutron 
flux with temperature and the delayed neutron concentration with time. 
The first assumption states that the reactor temperature rises instantan- 
eously with the flux (this is the usual prompt feedback model). Equation 


(7) becomes 

AE aa 
ar op 
aR _ a 
OT a Gee 


The second assumption states that the delayed neutron concentration 


and 


P 


has delayed neutron concentration has reached its steady-state value. 


15 





TABLE 2 [1] 


Delayed Neutron 


Group Half life (sec) 
Bro! 94.5 

Tee 24.4 

a 02 

re $ 

Bre 4.4 
Rp93: 34 e 


The justification of this assumption lies in the comparison of the neutron 
life time in Table 1 with the precursor half-life in Table 2. It is evident 
that a number of neutron generations will nave progressed setore the pre- 
cursor concentration has been changed. 

The assumption of steady-state precursor concentration reduces 
the neutronic equation to the case without delayed neutrons. This work, 
therefore, deals with the case where either no delayed neutrons exist 
in the reactor or the delayed neutron concentration instantaneously follows 
the flux distribution. The latter case is realistic when considering the 


reactor response to small perturbations. 


16 





IH. ANALYSIS OFTHE RCEACTORTWITH A 
TEM ERATURE COEFFICIENT 


— 


a re 
dT i 


(T”! Model) 

The solution of the nonlinear equations, Eqs. (3), (4) and (8), 
governing the reactor kinetics is accomplished by means of the Poincaré- 
Lindstedt perturbation technique [5]. This method offers a solution to 
the neutron flux distribution for small perturbations in the flux intro- 
duced by the temperature feedback effect. The perturbed flux distribution, 
referred to as P(x, E) is described by a Taylor series expansion 
about the initial state P(x) where € isthe amount of perturbation 
from P, Cx) v Te perturbation technique yields multiple equilib- 
rium solutions for Px E) and associated eigenvalues for each 
of these states. The presence of more than one equilibrium state, known 
as the bifurcation phenomenon, requires the investigation of the stability 
of each of the possible states of flux. Stability considerations are made 
by imposing an initial deviation on an equilibrium state and examining 


the resulting temporal behavior of this state. 


A. EOUTLIBRIUM STATES FOR ATSLAB REACTOR 
Before the general solution to Eqs. (3), (4) and (8) is undertaken 
the equilibrium states for the slab reactor is investigated. For the case 


oh 


where the flux does not vary with time (i.e. DE = O ) and the 


17 





reactor kinetics are not dependent upon delayed neutrons the governing 
equations reduce to the following boundary-value problem in the spacial 
domain D[-1, 1): 


Vb + (bly) -1)9) = 0 mo 
Pix) = O a SD 


and 





— 


A R(P) _ O 
ap e 


Rip) = RR: 


Equation (8) may be solved to give 


Pox) 
R (e) = R - a Êm TES (9) 


This result may be substituted into the boundary-value problem and the 
result is 
2 C(x) 
V W(x) A (> =O Im TREF | 00 = 3 in 


(10) 
Qx) = O on SD 


where N = R. = { 





U 


Tie solution to o os constructed about a flux distribution 
P, (x) and eigenvalue A o which will be determined. 
Physically, A o is this eigenvalue corresponding to the initial 


material buckling of the reactor before perturbation by feedback. After 


18 





the system is perturbed the flux and the material buckling establish a 
new state which may be described by a Taylor series expansion about the 
initial state. If the amount of perturbation is measured by € , the 


resultant reactor state is 


P(x) = Px) + € Pox) + LE PIX) E 


Co 


E A A a a 


(11) 


Thus for small perturbations of the system the values of Pix, €) 
and ACE) may adequately be described by truncating the series aíter 
the third term. The unperturbed flux, op, (x) ‚ and eigenvalues, 


A ES , are defined as the solution to the linear problem 


IP + D a) Y (x) = O D 
Rx) = © on 6 


Since Ve and Q are constants,the solution is 


Ctx) = Alm Cos YA a xX 
2 
k-a = ee 


where Mr 1, 3, 5, ..- . For the purpose of simplifying the subsequent 


and (13) 


analysis the shape function W for the slab reactor is denmicd a2 
MI ) 
v, = Cos (" x (14) 


19 





The remaining terms of the expansion given by Eq. (11) are determined 
by differentiating Eq. (10) with respect to € and taking the limit as 


€ tends to zero. The values of p and A are then determined. 


By repeating this process the values of op and A are evaluated. 


Maas, differentiating Eq. (10) once yields 
2 A o y J Y e 
FT PAP + AP- OA ln YP- aY =o, (15) 


The limit as E tends to zerò requires that 


fim Pexe) = P ex) 
€ > O 


Kim ACE) = N $ 


E -> 0 


As € tends to Zero, La. ISl redu cec to 


T Lex) ze Ea a) WES = ALO in D (16) 


(D(x) = © on SD, 


If Eq. (16) is to have a nontrivial solution, the solvability condition 


(provided in Appendix A) may be evoked. Thus 


| -À Poo) Pix dx = O 
D 


un 
where P (x) is the homogeneous solution to Eq. (16). The 


' H SN H al 
solution for co (x) 18 easily seen to be Y (x) F Br ce i 
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amily, 
When the functions cp (x) and E (x) are substituted 


into the solvability condition the result is 


i H 2 
E D 
Since the integral is nonzero, it follows that A =© . This result 
is substituted into Eq. (16) and the solution for Op (x) is clearly 


seen to be (x) = Bn A . Thus the linear terms of the expansions 


are 


Y lx) i a vy (X) and 


= © 


The parameter E is chosen to represent the average magnitude of 
Mie Jinear term in the expansion of Eq. (11), weighted by the mM 


distribution, and measured in units of the maximum initial flux. Thus, 


l 
€ = An | Kew ( (Pix, €) - 2.0) dx . (un 


Bifterentiating Eq. (17) with respect to C yields 


An = [on Pix e) dx 
D 
where 


Pexe) = Pixi + E Pix) + 


ZA 





By letting € approach zero 
2 
An = Bul TY ~ dx 
D 
thus Am = Ba . The second terms of the expansions are 


Y = ae Y, and A = JO. (18) 


The third sem of Pix, é) and AC&) © are similarly 
determined by differentiating Eq. (16) with respect to É . The result 
is 

Pine) + AW Pine) + DI PHO 
+ Ale) Pine - a. Pines In (TE) 


PCR) 
- O PIO Pra — a Pxée) = © in D 
P(x = O on SD, 


Fince A= O7 the resulit-of talong the limit asm c Ferenc torzero 


is 


Co ot 
V Poo + (A.- a) POX) = (19) 
. 2 | ./ i 
A Pox) Pix) 7 A P cx) in D 
24 
tx) = O on SD, 


The solvability condition requires that 


e Z | a ee H 
lack (Dox) 7 Na Oxy dx =O 
D 
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fe H 
where p (x) is the homogeneous solution to Eq. (19). This 
sı H _ h p 
solution is obviously seen to be P (x) 7 E mM 
When the various functions are substituted into the solvability condition 


the following becomes evident 


A An ch] $, dx = a An C J Y da 


and A OE . This result is replaced in Eq. (19) and it is found 
that the right hand side identically becomes zero. Therefore, the solution 
for the third terms of the expansions is p On a= Ge wv, and (20) 
\ =< (O... The cgquilibrium solution for Px, E) is obtained 


pa placing the results of Eqs. (13), (18), and (20) into the series of 


eee iil). Thus, 


Pex = Ay (146 +80 Ste) oy en 


and 


2 
K-a = (“,) +7 eat 


Mee MM = 1, 3, 5... + The value ot Pix €) from Eq. 21) may 
be evaluated if the constants A, and C} are defined. Since the reactor 
is assumed to operate at its fundamental mode before temperature feed- 


back becomes effective (later it will be shown that the other modes are 


unstable), the value of A, may be determined by the power requirements 


23 





of the reactor. Thus, 


P = | (Cx) ax 


where P is the linear reactor power level and E is a conversion 
factor. The value aña is undetermined and will remain an unknown 
2 

of order € En throughout the analysis. For each mode, ™ , P(x, €) 
Pemeeciined in Eq. (21), describes two possible states in the reactor 
as depicted in Figure 2 for the fundamental mode. The state of larger 
flux is determined by a positive value of € while the smaller state 
is described by a negative value of E . Associated with each state 
of the reactor is an eigenvalue, ACE) = CS "given ber gar) and 
plotted against € inFigure 3. The values of AC€)- A forthe 
case of negative reactivity are drawn as a solid parabola which opens 
to the right. It may then be noted that the equilibrium state may only 
exist for values of A (€) =" Come which are oreater thanithe veque of 

A o «+ the initial material buckling of the unperturbed system. For 
the reactor with positive reactivity the converse iS apparent. The cor- 
responding curve for positive reactivity in Figure 3 is that drawn with 
a dashed line. These eigenvalues are found to be always less than the 
value of Ao . It may also be noted that the eigenvalues of the 
System are ALE) - Q where the ao — A  ) describes the 
amount the material buckling changes so as to accommodate the temper- 
ature feedback. 
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PLUX DISTRIBUTION FOR SLAB REACTOR 
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B. STABILITY ANALYSIS OF THE SLAB REACTOR 

The stability criterion for a slab reactor is developed by considering 
the temporal effect of an initial deviation on p x E) . Under 
these conditions the governing equations, Eqs. (3), (4) and (9) reduce 


to 


RI E A: 


Pixir) | 
“= > in D 

Ro) EF Ro Qu. Ln (x) (9) 

Y (x,0) une) j P(x, t)= O on SD, 
These equations may be compacted into the more convenient form 

OY Cx, T) Poo 

A = Y Px, 7) +j] A- O Un e ay V(x, 7) 

in D (22) 
CP X) T) = O on SD 
where 
A= R- 1. 
The time dependent nature of the neutron flux is investigated 

by prescribing an arbitrary initial deviation, H , on the equilibrium 


flux and then observing the growth of the.deviation. An unstable state 
of flux is that which will result in an unbounded growth of the deviation. 


A stable system will behave in a neutrally stable manner which implies 
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that the deviation tends to a new equilibrium state or in an asymptotically 
Stable manner which implies that the deviation decays to zero. The time 
dependent flux is described by a Taylor series expansion about the 


equilibrium state (x, E) , defined by Eq. (21), 


(23) 


} 2 
Pine) = Pue + N PRO 42 Pog GED 40 
such that 


Rim Pix € tn) = 


->o 
Equation (22) is differentiated with respect to EC, Therefore, 
d ia = 7e MNAE 0 raS ) Y, in D 
DT 2 E ~ ESOO A A 
a = O on SD 
where 


h = oe (Ke Ty) , 


The limit as u tends to zero is taken and the result is 


Oo Fy, CX, €,7) 2 
re = V E, (55) (Ace zu (24) 


(x, 7) ; 
al Un Ze) E, Ox ED) In D 


A (x € 7) = O SN 
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If a solution for a (x, És T) is chosen to be 


m = UKE) exp ( olce) T) (25) 


then the stability of the system is analyzed by inquiring whether the 
valueof &X(€E) is negative or positive. A form for ol (€) 


is given by the series 


of) = Kr ed dem: es 


The value of A (€) is determined by the perturbation tech- 
Mique previously used. In this vein, Eq. (25) is substituted into Eq. (24) 


and the resulting equation is 





2 
NY U C(x E) + (ce Ml CES = (27) 
P(x ,é) | 
Ou Im pie) | LCi) = ©) in D 
uU (xE) = O on SD, 


Setting E equal to zero 
2 (28) 
Y U(x E) - (Aon e ol.) Ulx,o) =O in D 


U(xXO) = O on OD, 


Since the values of A om ecu Q. ) have been prescribed by Eq. (13), 


29 





Eq. (28) may have a nontrivial solution only if (proof is given in 


Appendix C) 
A om Q = A ee CR ol o 
Thus the value of Oo, is defined as 


K Be (Xom- a) - (Aom- a) (29) 
= (Y) Um? - m) 


for odd intergers of M and Mi where ym is the mode of the 
deviation. If the fundamental mode of flux distribution, Mm = E i 


considered, it is evident that Eq. (29) reduces to 


oim = (1) (1- më) 


ime the values of Xo, FR are 


Kom = O Ma = He 


30 





For values of € near zero the first term of the series for CE) 
indicates that the fundamental mode is not unstable for any value fm 
However, the neutral stability predicted by yen =] obliges an investi- 
gation into the higher order terms of Eq. (26). Some conclusive results 
may be obtained by considering the linear term of Of (€) 


Hor highenmmodes o y = 1 ) Eq. (29) indicates 


that 
ol (€) = = 0 mn 7 mM 
= © rn =F IW 
= 2 
= L mn “PR 


Since a positive value of ol (€) may OCCür, it is Concluded that the 
higher modes of the reactor are unstable. For a physically realistic 
Situation, however, only the fundamental mode exists in the reactor 
before the perturbation. 

The solution for U CS O) is obtained by substituting the 


Wale of Eq. (29) into Eq. (28). The result is that 


2. ; 
7 uno) + (Aom- A) UK) =O) in o 60) 


RT A 


lt is then obvious that 


OS e 
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where the index mm designates the modes of the deviation. The value 
[N= 1 refers to the fundamental mode of the deviation while [in > 1 
refers to the higher modes. 

The linear term of the expansion for œ C€)  , defined in 
Eq. (26), is derived by differentiating Eq. (27) with respect to € 


thus 


lo > 
Y IES e (Ace) - a - ole) (31) 


Pix En. ' ; 
on od owe + (her - dee 


W(x, €) 
- Q Eh U CXE) == O in D 
CEO n SDi; 


The limitas € tends to zero is taken and 


(32) 
TV" UCK,o) + Be E eae) U (xio) = 


| l Pox) | 
E ACo) + o(o + @ ad LACX,O) in D 


ú (xo) = O TET 


By noting that x =n) the solvability condition for Eq. (32) 


requires that 


. x ; 
| + O 5, ) U (xo) Lo) dx =0 
D O 


. H 
where LA ( x,0) is the homogeneous solution to Eq. (32). When 
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the value of Ces , from Eq. (29), is employed the function 


e H 
wy: CX,0) is determined to be 


e H > H 


The value of © is defined by substituting the functions P (x) ; 
a o H 
P (x) (rom Eqs a anda ELA) , and ATX 0) 


into the solvability condition and the result is 


AR Olen. (33) 


This result and Eq. (29) determine the linear approximation to 


SUE) , Eq. (26). Thus, 
e 
HCE) = (Tr) (m -mi) - Ea 


It now becomes evident that the question posed concerning the stability 
of the fundamental mode ( [Ve ) for sufficiently small € is 
answered by the linear approximation for œ (€) . Since values of 
MM > Y give very large, negative values of A (€) and thereby 
assure the stability of these modes, the concern is centered about 


m = mM = 1 . For this case 


AUS = — Ea 


and the stability of the fundamental mode will depend upon the sign 
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of the linear term of oA ( €) . The various cases to be analyzed are 
the following: 


a) negative reactivity 
of ( E) SHO) for positive € 
> O for negative € 


These equilibrium states of flux are depicted in Figure 2 and the corres- 
ponding eigenvalues are shown by the solid curve in Figure 3. A positive 
value of E describes an equilibrium state of flux which is larger 
Maine initial flux. This state, Eq. (21), exists for an eigenvalue, 

A le) —- Q__, greater than the initial material buckling of the unper- 
turbed reactor, Ào , and it is observed to be stable. A negative 
value of € indicates an unstable equilibrium state in the reactor. 
However, Eq. (21) indicates that the eigenvalues must be greater than 
the value of Au for any value of q . Therefore, a situation in 
the reactor which permits a decrease in flux (i.e., negative & |) with 
an associated increase of the reactor buckling is physically unacceptable. 
b) positive reactivity (evaluated by changing the sign of all terms con- 


tenio CL) 
xL (€) > O for positive € 


< O for negative € 


34 





These states of flux are depicted in Figure 2 and eigenvalues shown by 
the dashed line in Figure 3 are analyzed in a similar manner. The stable 
state of flux is that described by a negative value of € . From Bee 
noting that the sign of Q must be changed to accommodate positive re- 
activity, the fundamental eigenvalue is less than that of the initial 
Material buckling, he  Inepei@re, 9a decrease in tithe (i.e. negative 
€  ) is compatible with 'a decrease in the eigenvalues. However, for 

the case of positive € the equilibrium state is unstable. Again it 
is noted that the increase of flux associated with a decrease in the 
material buckling is uncompatible and this situation is physically un- 
realistic. 

The solution to Eq. (32) can be completed by first substituting 
Eq. (33) into Eq. (32). The right hand side of this equation becomes: 

Dl xX) 

RHS = Cto) +a xy) Ulx,O) . 

However, from Eqs. (13), (18), and (30); P ix) = A m WV, > 


P(x) = Ave, W, ) and UL(x,0) = Un respectively. Thus, 


Rus = (-a + ac Un Yo 
= Sr 


Equation (32) correspondingly reduces to the form 


39 





Two + (Aom- a - hom ) UX,0) = O ind 


Ulx,o) = O on $D. 


Thus, the solution to OXO) is 


ow : 
U (XO) = ex) = Us Wn 5 (34) 


For larger arguments of E 


(E) 


the third term of the expansion of 


, Eq. (26), may become important. 


This term may be evalu- 
Bee by differentiating Eq. (31) with respect to E 


and then taking the 
limit as E 


tends to zero. The differentiation yields 


TUE) + (Ace = 0 - &(6) 
a An ee) U (0% E) +2 (Ace 


= AO 0e 
OGG) hae Pa) LETS) 


+ (‘le = lee dee AS 





. 2 ae) 
(0 nt. E 

+ 0 pri) ute =0 n a 
U = © on éD., 


When the limit as E tends to zero is taken and the various known 


functions are substituted into the resulting equation, the following 


becomes evident. 
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TI” ik Gx,0) + rom — A- Komm) Ú (x0) = 63) 
lk - a(2- 5 )| Um Ym in D 


Lt (%, 0) = O on SD. 


The solvability condition requires that 


[ [& - 0. (2-FE)] Un Ya ittsodx <0 


where \A (%,0) _ is the homogeneous solution to the above equation. 
oe H Ar H > H 
The nontrivial solution for UA (X,0) is Lk €x,0) = O WPa 


Thus, 


ER a (2- Af a, tim Ape ANO 


2 
and since the ae A dX = 1 ’ 


de re 
= eee) 36 
ol a(2> 36) 
The expansion of ol le) , Eq. (26), may be written in the more 


evn cte form by the substitutions of Eqs. (29), (33), and (36). Thus 


%) 
A (cd) =(% (m? -m) -€4 
Grn 
+5 € A (2- Rote (37) 
For larger arguments of E the stability of the system is contingent 


upon the negativity of ol (€) Maetinea Dv Ed. 3/). For the 


Sí 





fundamental mode ( M= Í ) Eq. (37) reduces to 
O =(%) (12 ml A 
ote) = 2. (4- m') -€a+7€a(2- A 
The value of fM of particular interest is that given by AM = 1 since 
values of M > li provide a large, negative first term in Eq. (38). 


Thus, the stability criterion is 
A 
ca + 2€0 (2- Ya.) < O 
Or 


ca (1-¢€+46¢€%)<o (39 


By inspection of Eq. (39) the quantity within the parenthesis may be 
seen to always be positive for small arguments of E > Thus, the 
stability is governed by the coefficient (- É A )as has previously 
been considered. However, for cases of larger arguments of € the 
inequality shown above may be changed. Thus, to the order of E 


the condition for stability is 


1 
E - Ib = z YA, | a 


If this were violated, such large values of E would render the re- 
actor unstable. Thus, the undetermined peak amplitude of Pix) j 
C} , aS compared to the peak amplitude of the unperturbed flux, Ay. 


is important when considering the stability regime. 
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oF EQUILIBRIUM STATES FOR A FINITE CYLINDRICAL REACTOR 

The finite cylinder represents a realistic geometric configuration 
for a reactor. The physical and mathematical considerations which have 
been deduced for the slab reactor are directly applicable to the cylindrical 
Heeetor. In this section, the equilibrium states of flux, associated 
eigenvalues, and the stability of these states which have undergone 
feedback perturbation from some initial state will be studied for cylind- 
rical geometry with azimuthal symmetry. The flux distribution must 
depend upon two spacial coordinates identified by 8 ,ı the radial 
coordinate, and £2 ‚the longitudinal coordinate. 

The equilibrium solution for the reactor is developed from the govern- 

O P Ce eae) ar 

ing equations under the conditions that y =O A 
The governing equations reduce to the form given by Eq. (10) which is 


rewritten below, 


(41) 


2 Pir, z2) 
JI Pinen +(a - A Am en) Pan DOR 


- P(r2,7T} = O on $D. 


The feedback-perturbed solution to Eq. (41) is evaluated by a 
Taylor series expansion about the initial state. If the amount of per- 
turbation is E , which will subsequently be defined, the resultant 


equilibrium state is 


39 





5 zur 
Wr, 2,6) = P (12) + E Praz) + E Pra) +... (42) 


Lu 


ACE) = AS + CA as zé A teese 


The initial state, £ Cr, 2) and As ‚is defined as the solution 


to 


= 
F Euz) + (Aa) Laza = O inp 


ROZ = O on SD. 
Equation (43) is solved by assuming a separable solution for L (G 2) 


given by 


Ar) = Rin £w. (44) 
Substituting Eq. (44) into Eq. (43) yields the following: 


2 
N Le u 
Rtr) a 2 T £ E) ae Av 


t (Ao- a.) Ai == O 
which may be separated into the equations 


ARCr) 2 
a fr torts Rir) = O in Dd 


dr Ar 
Ror) le on óD 9 
da Az) 
“3 + RA 0 la Dd 
FI) = O on ë óD 5 
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The solutions to these equations are the following: 


Rn = Am dp (an a) > 
Zar = An (los ("M 2) ; 
Mi >.-q = (My + ("IA y 


where Arm are inet or mmerzeroin order Bessel function: The 


Solution to Eq. (43) is 


Powe) = Rn) Zoe) 
Am Je (um 7) Cos ("I z ) 
mm ra = (Moe) ty 


N 


h a re 
where ME 1 ye 


For the purpose of simplifying the subsequent analysis the shape function 


Sa is defined for the cylinder as 


14 = J, (Ma Zr.) Cos (MI z). ao 


Thus, 


Pa eee An ap (47) 
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The remaining terms of the expansion of P ( (ace €) and 
A (€) , given by Eq. (42), are determined by differentiating Eq. (41) 
with respect to É and then taking the limit as É tends to zero. 


These steps lead to the equation 


A . . 
Y Perz) + làa) Pinza =- AP np 


Porz = O on Sp. (48) 


If Pir t) is to have a nontrivial solution, the solvability con- 


dition (provided in Appendix A) may be envoked. Thus, 


: AH 
e AL (12) Pie 27 r drar+ =0 5 
D 


where P = 2) is the homogeneous solution to Eq. (48). This 


solution has a similar shape function as Eq. (43). Thus, 


-y H 
P (Pa? — B Vy a (50) 


ie 
The fundtions P or >) and P Gp Has. (47 ioeang (50), 


respectively, are substituted into Eq. (49); hence, 


f -A An Yo Bi Y, 2y Y ana Zz. — © 


- A Am Bn 2y) W fdr dz = 0. 
D 
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Since the integral 
er e 
J y r dr dz = 2 ie o (Mm) 
EOR 


the solvability condition may be satisfied only if A =O . This 


result may be substituted into Eq. (48) and the value P (r, 2) is 


the solution to 


J‘ Pira A, a) Pina) =0 in D 
Pine) = © on SD. 


Hence, 
3 E 
7 


Pine) = Pines = Bm Y. (60 


The parameter € is chosen to represent the average magnitude 


of the linear term in the expansion of Eq. (42), weighted by the function 


_ Vn _ 
yi lo J Gam) 


Thus, 


and measured in units of maximum initial flux. 


| ab | (51) ä 
_ m (6 2€)- a)? rar dat 
E Am) Tr? Syn) + ta] 





Differentiating Eq. (51) with respect to E and letting E approach 


zero yields 
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i 


Bn] Me: 2yr_ dr dz 


An er I (Hm) 


= Bm. 


Therefore the linear terms of the expansions of Y ( FEP E) and 


NCE) are 


P = Am Y, 
and NO 


The third terms of the expansion of Or Č, E) and N (€) 
are determined by differentiating Eq. (41) with respect to E twice 


and letting É approach zero 


Y Pira + (A-a) Pira = (54) 


Pina y Pier) eo.» 
C 
CP (r, 2) = © on SD. 


The solvability condition requires that 


Bine y 
(a GE, ~X Ree) Poon apr are =o 


where Y ie 2) is the homogeneous solution to Eq. (54). The 
-> H 
solution for Y Cr 2) is similar to that obtained rom Eq. (43), 


thus 
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. H H 
Pinza = Cy ur. 


The functions Y Cr 2) and P (1, 4) , determined from Eqs. 
= H 
Cand (53) respectively, and Y Cr 2) are substituted into the 


solvability condition. Thus 


| (An Im - Dd An Yn) Cn Vs, 2a de p 
5 (a-3) [| Per de dz <0. 


Mreretore, 


A iS QA. (55a) 


If this result is placed into Eq. (54), the right hand side will identically 


be equal to zero. Thus the solution for p Ge / z) is similar to that 


Obtained for Eq. (43) and 


Pr) = CN we (55b) 


The various results for the terms of the expansion of Eq. (42) are 


presented as follows: 
Pleme =A, (146 +56 e+) mH, 
ao- a (Pr) (MY) + ¿Ear 
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It is apparent that for any mode of flux there will be two possible 
equilibrium states defined by the negative and positive values of é 
For each of these states a unique eigenvalue is present. For negative 
reactivity ( Ale) - AL) must be greater than Ao and for positive re- 
activity the opposite must be apparent. It may be noted that the eigen- 
values are defined by ACE) — Q . Where the tem -A denotes 
the amount of change of the material buckling to accommodate the tem- 


perature feedback. 


19% DTAP ANA NoE OF TAHE FINITE CYLINDRICAL REACTOR 
The stability criterion for a cylindrical reactor is developed by 
considering the temporal effect of a deviation on the equilibrium states. 


The governing equation has the form 


O Pon ns 2 Pirat) (57) 
E = VW Pirar) +()- A Am a) 


i" D 





Por, rt) = O on 6D 
A CE 


The temporal nature of the flux is investigated hy imposing an 
initial deviation, N , on the equilibrium states defined by Eq. (56) 
and then observing the growth of the deviation. An unbounded increase 
will indicate instability; neutral stability implies that the deviation 


tends to a new equilibrium state; and asymptotic stability is assured 





by the decay of the deviation to zero. The deviation is described by a 
Taylor series expansion about the equilibrium state Per £, E) , de- 


fined by Eq. (56): 


Pr 2, EN, T) = Pinze) + N A (GET) (58) 
\ è E 
+ zn PbO BET) = = 


such that 
fim Pirae) = mao 
YO 


Equation (57) is differentiated with respect to q and the limit as 


u tends to zero is taken. The result of this process is 


(58a) 
ò ta _ NN) Y (AW -aA - a Im EBON Ooi m D 


% = Pa (2E) A ZC) on Jd 


If a solution for pn 2,€,7) o e 
, (12,60 = (U(NZ,g) exp ( oe) T) 3  (58b) 


then the stability of the system is determined by the sign of Ale) 


OLI LOL al (€) is given by the series expansion 


AU) = & + €& ute ol +... (59) 


The value of of (€) is determined by the perturbation technique 


previously used. The function given by Eq. (58b) is substituted into 
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Eq. (58a) and the resulting relation is 60) 
Per az,e ?, EN 
=0 


Vu + (ao - Blei a im Der) 


im D 
(A(G 2,6) = O on JSD 
The value of Xe may be determined by taking the limit of Eq. (60) 
as E tends to zero. Thus, 
2 (61) 
Yume + (A, Aa -0,)uan»=0 ind 
on dD, 


Ulnz) = O 


A nontrivial solution to this equation is obtained only if (see Appendix 


Caor proof of this) 


AS SN = peek - QA- lo : 


Thus As is given by 


Oem ae e A) = A ne a) (62) 
A Uan pum) (AY (mm), 


Some remarks may be directed toward the stability of the system via 


li 


Eq. (62). As was apparent for the slab reactor the fundamental mode 


( m =] ) may be evaluated as follows: 


Ones = O MM = 1 (62a) 





Thus the higher modes of deviation (met ) which are imposed upon 

the equilibrium flux will vanish with time. However, the value oí m = 1 
provides A Ó and the system may be thought of as neutrally 
stable. Thus it is necessary to inquire about the effect of A upon 

the system for the value oí [IM = 1 ? 


The higher modes of the reactor [ M> 4 ) are analyzed as follows: 


xK (E) < ~25y" NO > Mm 
= O m em (62b) 


2o vn G MNA, 


Since the condition provided by mM «= N may be realized for modes 
of M _ > A , values of x &E)> O are present and the higher modes 
are unstable. However, the reactor is assumed to be operating at the 

fundamental mode before the feedback is introduced. Equation (62) for 


ol omm may be used in Eq. (61) to provide a result similar to Eq. (43). 
Thus, 
wird = Ur Yon. (63) 


Bee Appendix C {or the prooi ol iis.) 


The linear term of the expansion of ol (€) (Eq. (59)) is evaluated 


by differentiating Eq. (60) with respect to €& and taking the limit as 


49 





€ approaches zero. The result is 


an ‘ 
N =) + E EA AS U (vr 2) (64) 


x Š (r2) 
n CO A pan) ULF Z) nD 


Poor, 2) 
Ber =l O on dD 


e 


By noting that A has previously been determined to be zero the 


solvability condition requires that 
lo + Pina A U Ling = O (64b) 
Ara Pira) (G UTH = 
D 


- Y 
where UL (r, +) is the homogeneous solution to Eq. (64). By the 


Same reasoning used for the solution of A(G Z) the solution for 
» H 


au! ae 


= 


When this result and those provided by Eqs. (47), (53), and (63) are 


used to evaluate Eq. (645) a value for A is derived. Thus, 


RE = A) O qa ¿Tr ar az=0 


(ro) Y ur q =O, 
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Therefore, 


a ae (65) 


This result may be substituted into the expansion for OA (€) and 


further analysis of the stability is made. Thus, 


œX (€) = ae == EQ A (66) 


The question pertaining to the stability of the fundamental mode (m = 1) 
is resolved for small € . For the case where ONE 1 the value of 


of oll is zero and 


Al) = -Ea . 


The analysis of the various cases of negative and positive feedback is 
similar to that developed for the slab reactor. The conditions of stability 
for the fundamental mode are briefly stated as follows: 

a) The stable equilibrium state for a reactor with negative reactivity 

has a positive value of E . The other en state is unrealistic. 
b) The stable equilibrium state for a reactor with positive reactivity has 

a negative value of € 2 (me sOtner Cquibibrimi state is Unrealistic, 

If the value of A is substituted into Eq. (64), the right hand side 
will become zero. Therefore the solution for U Cr z) iS thE 


homogeneous solution to Eq. (64) and 


iL (4, 2%) = Um N m > (66a) 
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For larger values of É the quadratic term of © (€) may be 
necessary. To obtain the desired results Eq. (60) is differentiated twice 
with respect to E and the limitas É  tends to zero is evaluated. 


If all of the known functions are employed the result is 
2 0 or 

Y ts sid Na: ia PA 67) 

Y Ce 2) = O On SD 


The solvability condition u that 


Jl - a(2- a ) | te Wa Ui ira) dar dy dz =O 


H 


ee 
where UL CG E) is the homogeneous solution to Eq. (67). Therefore 
ta Fr. úl H NV 
. It is obvious that the solvability criterion is met 
only if 
© v om 
or 


| 


o A A (68) 


This result is substituted into the expansion for ol (e ) ES (59), 


and the following becomes apparent, 


OZ 


69) 


2 ( 
RE) = I) (un 53 HAY (mi má) 
-ca +7 EA h2- VAn) te 


Since the fundamental mode of operation is generally of greatest interest, 


Eq. (69) is rewritten for /M = dle . 


x E) = T a Me +) (1- m) vo 
z éal 4 -7€(2 - “Ya,)| fore 


For PM > 1 the first two terms dominate and A (€) < O STO 


pm = 1 the stability may be determined by evaluating 
L C 
xE) = -a| 1- ze (2- YAn) - 


These results are similar to those obtained for the slab reactor and 


Stability is granted under the conditions established for Ea. (37). 
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III. ANALYSIS OF THE REACTOR WITH A 
TENIFERAFURE COEFFICIENT 


Ze 302 


ER (y MODEL) 


aT 


The solution of the nonlinear equations (3), (4), and (8) governing 
the reactor kinetics is obtained by the same perturbation technique used 
to analyze the pl model." This method investigates the effect of a 
Small perturbation introduced by the feedback effect. The equilibrium 
state that ensues from the perturbation is defined by a series expansion 
about the initial unperturbed state of the reactor. The perturbation 
technique offers multiple solutions and unique eigenvalues for each of 
the equilibrium states (see Figure 2). The presence of more than one 
equilibrium state, known as the bifurcation phenomenon, requires an 
investigation of the stability of each of these solutions. The procedure 
of analyzing the stability entails imposing an initial deviation on the 
equilibrium flux and then examining the resulting temporal behavior of 


the deviation. 


A. EOUITLIBRIUM STATES IFOR TAE SEAP REACTOR 

' ‘ e ‘ + E] -3/2 
The first consideration to be made in the analysis of the "T 
model" is that of determining the possible equilibrium states of the 
reactor after a feedback perturbation has been imposed. The solution 
to the governing Eqs. (3), (4), and (8) is desired for the situation where 


OFT) _ 
the perturbed flux has achieved the equilibrium state (i.e., oK =0) : 
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The governing equations reduce to the form 


T Pew) te (P(p)-1) P(x) = O in D (3) 
P(x) O on AD 


and O 
dẹ B 


= NA E p T= 
ol op p /z ) R CP.) (8) 


The equation for the multiplication factor is easily shown to be 
a =e: 

Ri) = A al | Pox - Px]. 
The function $ (x) and the constant fe. are cOnsidered the 
initial flux distribution and multiplication factor before the feedback 
effect becomes significant. The result is substituted into Eq. (3), 
yielding thereby: 

) 
2 aie me (71) 
Ton +L ara ( Pex) - Px) (x)=0 ind 
Pix) = O on «SD 


where 
After the perturbation due to feedback the new equilibrium state 
is described by a Taylor series expansion about the initial flux distri- 


bution, L Gx) , and the initial eigenvalues, AS . If the amount 


of perturbation is measured by € (the perturbation parameter which 
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will subsequently be defined), the resultant reactor state is 


5 basen See 
Pine = Pod +€ Pon +7 Pia +: ay 


A p Laa 
IAS E TEA FEA 
For small perturbations of the system the values of Plx,e) and A (€) 


are described by truncating the series after the third term. The un- 


perturbed states for the reactor is defined as the solution to the problem 


È 
V x) + 2 Cn) E in D (72) 
Pix) = 0 on dD, 


It is readily seen that 


= A m Cos VA.” * 
kl CU (721) 


m = j 5, 50 ng a) 
So as to simplify the subsequent analysis the shape function m is 


defined for the slab reactor as 
Vn, = Cos (“T/ ) De va 


Thus, 


LC = AR A f (74) 


The next term of the expansion is determined by differentiating 
Eda. (7) with respecto E and tooling che limitas E tends to 


zero. This differentiation yields 
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: | yh _ (75) 
T Lixo t | AC) au a (Pexe - P cxe) P(x, €) 
OUS 2 
ze ( NE) 2 Poe Pexe) Pixeo in D 


Pixe) = O on dD. 


As é approaches zero, Powe) approaches DP ox) : A (E) 


approaches Ao ‚and 
a Zur Klee; (76) 
V PX) + O. es = P. F) Pix) = = A Polo) 
p m D 
Pix) = O on dD, 


inerder to obtain a Solution to Ea. (76) the solvability condition (given 


in Appendix A) requires that the 


: . H 
| Cig) For, Wax =O (77) 
D 


where Box) is the homogeneous solution to Eq. (76) given below 
yo Dix, T DN y 5 R a BL a =O me 
. H 
P (xi = O0 on SD, 
A solution to Eq. (78) is offered by a modification to the R. Courant 
and D. Hilbert perturbation method (provided in Appendix B). Their 
approximated solution is given by 


> H G oe 
Pix; = LL Cx) fe E ex + (2 ) Er 


a 





2. 
A 
Since the contribution ol ( /2 is very small only the first two terms 
cmne series will be retained Ihe substitution of Eq. (79) into Eq. (78) 


leads to the requirement that Umoo be the solution to 
T 
Umtx) = O . on dD 


the solution to Nm (x) can be obtained as 








09 
aie; 
Nmtx) D A mes Uy 
j om Y 


Jz USS: 
J EM A 
where dmi = P ix) ER A, ol x 
The orthonormal solution for UmnlX) is Unlı) = WPa and 
& 
Dem = (IA ) . The constant Al my is 


determined as 


| R 
Am; 7], Y Vy AK ; 


09 ft of E | 
UCA) - E m? i Y, 


J= 1,3,5,: ey 
"m 
These results may be combined into the expansion given by Eq. (77) to 


thus 


yield the following result, 


en H 2a e Glen 
Pix) = Br K a J me- gE kr pa (80) 
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For the convenience of the subsequent analysis this function will be 


referred to as the cig function given by 








= W + A a , (61) 
JE ATS J l 


A question may arise di the convergence of the series 
used in Eq. (81). If the Weierstrass M test for uniform convergence is 
employed, a series which is a continuous function of X inthe domain 


D is convergent under the condition 
f | 


where y M is a known convergent series of positive con- 





dants [6]. If the convergence of 


Fy (x) = dmi 


m= 4 


| 


N 


Y. 
Ty tx) 


is doubted this test may be used. First, the maximum value of 


is determined as follows: 
np? 
Y. AX 
A my y A, -f 4 e 
= 7 
dl % MAX AA AX 


e 
Um. 


IN 


IN 
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aa 
Since the function ae = 1 throughout the domain, 


LAA An, 
VA, (j°-m) ps mè- j? Y; (x) 


If the convergence of the left hand side of the inequality is provided, 
then f; (X) is convergent by virtue of the Weierstrass M test. The 
integral test provides a convenient method to establish this convergence. 
Hence, if the integral 
Fe 
| VA m (y m) dj 
o 


MAI 
exists, its convergence and that of J 2 ny (X) is determined. 
pas) J 





Thus, 
T EE 2 jz en 
J, An (em) Of = Jaz (7l Coth Gh) 
J 
a! 
> (+a E Cot) 
where E) is any possible value of j . Since the integral exists, 


the series converges. 


- H 
The values of p (x) and P (Xx) from Eqs. (74) and (80), 


respectively, are substituted into Eq. (77) and A is evaluated as 


shown. 
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JA An Y Ba | BY mein a ax =0 


A m 


Because of the orthogonality of the vv, function the above equation 


reduces to the following form, 
X H 2 . 
J-> An Ba Wi dx ="0. 


This equation may be satisfied by the condition 


A ZAC (82) 


The value of A is substituted into Eq. (76) and the homogeneous 


equation that results is similar to that defined by p Cana). Thus, 
Pex) = Ban Tr. (83) 


The value of the perturbation parameter, E , 1s chosen to represent 
the average magnitude of the linear term in the expansion for Eq. (11) 
weighted by the =. function and measured in units of the maximum 


iita! flux. Thus 


l 
é = zul, T x [ Pese) =P, cx)] dx (84) 


If Eq. (84) is differentiated with respect to E , the result is 


| e 2 
Ps An] Lo L Por + € Peor lax. 
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By letting € approach zero and employing the solution for Pox) 


Mes, Pix) = Bun WG ) the above equation reduces to 


Am 


2 
Bm sl de dx 


or, from equation (81) 


Am _ 2a Im; | 
Bm he + HD oe 





The quantity within the parenthesis is expanded and 


A 2 20. — Ama 
= = / an —— , 
B m Ed an + e] 7° ab af Ax 


m>- ji 4 
ee 
of 0° 2 
2i g DN 
+0 (5) me y Y * dx + O(a’), 


I TIS 
£ m 
The second term on the right hand side is identically zero because of 


the orthogonality of the Ab, function. When the integration of the 


a 
first and third terms is evaluated, i.e., d ap; ol x =i , the result 


is 





A 2/2 Am; i 
9=1,58,5),*"* 
Ir. 


E 
Since the order of QO. is generally very small, the contribution of 
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the higher order terms are neglected. Therefore, 
SNA 
A IM — B e 


Equation (83) that defines Pix) may be rewritten as 


Piw = Am T, (85) 


The third terms in the expansions for Pix é) and À (€) 
given in Eq. (11), are determined by differentiating Eq. (75) with respect 
to E and then taking the limit as € approaches zero. The dif- 


ferentiation yields 
Cre =e .> 
Y Dene + LACE + olPas - LEN) Poe 
A. 
+O or A Peixe P (x4) ) P (x €) = 
| a y 
E ( A (IA Pix e) P (x, €)) P (xe) 


- (Ye) - 


NID 


= 7% Sa > 72 
(2 (x,6) Poxey + 4 $ex e) P ixe) 


x Poxe),. 


As É approaches zero and Eq. (82) for Ato) is used the above 


equation simplifies to the form 
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T y 2 
T° Pix) + (r. j > Pon) & (x) = (86) 


-3/2 


es. Ena 
= A Pix) + S Px) P tx) in D 


Px) = O i on SD, 


Mae solvability condition for Eq. (86) is 


* LL 5.2 PH 
[C5 Rx) t S P ex) Pex) Pox) dx =O (87 
D 


au. è H 
where Cp (X) + the homogeneous solution to Eq. (86), may be seen 


to have precisely the same shape function as the solution to Eq. (78). 


Thus, 
Bi: Y 4 
Pox = Cy T a (88) 
¿ ee H 
If the values of Y Cx) i Px) , and Pe x) from Gas eur), 


(85), and (88), respectively, are substituted into the solvability con- 


dition of Eq. (87), the following becomes evident, 
“u H NY if 3 
MAC) Y, o AK = 4 FA, Er y AX, 


The integrals are simplified using the definition of Te given in 


Pa: (81 Thus 
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D 


because of the orthogonality of the Ab, function and 


[fee] ax E: y A A ax 


=1,3,S,** 
Em 


If the contribution of any term having a coefficient of as or larger is 
° 6 


considered negligible, an approximate value of A may be evaluated. 


Biene, 
< H RO N Ve 
A Amen = 4 Am Gm f Pa AT (89) 


The definite integral may be evaluated for the fundamental mode of 


è 


reactor operation (i.e., /¡W= 1 ) and the value of A for Me A is 


Q 
= A (.279) . (90) 


e. 


The equilibrium solution for Pix, éE) and A CE) is obtained from 
the results of Eqs. (12), (1/41, 182), (865). and (90). The expansion for 


Pix, e) and A CE) Mequation (las 
Pome = An (hte % + Oe) 
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where Mv, and MER are defined by Eqs. (73) and (81), respec- 


tively, and for the fundamental mode 


> Q 
ao =h + 7 ym G279) + 


2 


a 
(Y) + ¿0349 ¡2 Ec +... 02) 


1 


The magnitude of A, is determined by the "power condition" 


P = a, Pix) dax 


where P is the reactor power level and $ is a conversion factor. 
Equation (91) implies that there are two equilibrium states for each 
mode of operation corresponding to the negative and positive values of 
E . For the fundamental mode this condition may be depicted as 
shown in Figure 2. For reasonable values of € and (Q.) the con- 
tribution of ee is essentially that given by Ay, „he impact 
of the series of cosine terms within the ER function may be seen to 
be negligible if the peak amplitude € Á 13 MAX of the largest 


S@nmiributor is considered. 


r 


acc Ay 
€Aismax = € o =g Y, pa 


ANA | 
The value of a 13 = Ue was calculated and for a very 
\ 


a 
conservative value of E=.1 and ia an OOL 
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6 FOUL a 
ET = ci £ 1 zit 


€ | 


For each of the equilibrium states there are unique eigenvalues as shown 
in Figure 3. Fora reactor with a negative coefficient of reactivity the 
z 

fundamental eigenvalue is drawn as a parabola centered at ( I ) 
which opens to the right. Thus the equilibrium eigenvalues are always 
larger than those of the initial state. As has already been deduced for 

-] 
the "T ~ model" any decrease in the neutron flux (i.e., negative E ) 
cannot be compatible with an increase in the material buckling. There- 
fore negative values of E for negative reactivity are unrealistic. 
For the case of positive reactivity the eigenvalues react in an opposite 
way as depicted by the dashed line. Values of E which are negative 
represent possible states of flux, but positive values of É are un- 
realistic. 

| T. 

When the eigenvalues of the "T model" are compared with those 

11 -3/2 11 , s . 11 al 11 : 
of the "T model" on Figure 3 it is noted that the "T model" is 
affected more by equal perturbations on the system. This result is con- 


sistent with that predicted by A. Reichel [3] as can be shown on Figure 1. 
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-1 
There it is observed that the line for the "T model" is in very close 
proximity to Reichel's line of maximum Doppler effect. Hence, one may 
conclude that the Doppler effect has a more pronounced impact on the 


92 


Ej £ 
i model” reactor than on the "T model" reactor. 


B. Donny Ayolo OF THE SLAB REACTOR 
The stability criterion for the equilibrium states of a slab reactor 
with negligible concentrations of delayed neutrons is developed as 


follows. The governing equations, Eqs. (3), (4), and (8) are combined 


into the equation 


OY (x T) È =y, 
AS = V Pixa) +) - A (Pixa) Ca 


=Y2 ~) 
- P cx>) ] P(x) in Dd 
(x, T) O on éD 


A 


The stability of an equilibrium state is investigated by prescribing 


i 


an initial deviation, 7 , in the equilibrium flux and observing the 
growth of this deviation. If the effect of the deviation increases the 
flux with time, the system is said to be unstable. If the effect is 
reversed the system is stable. 

The time dependent flux is described by a Taylor series expansion 


about the equilibrium state, Px, E) , detined by Eq- Ol): 
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Deren) = Pie +” ye (XET) (94) 
| E 
ZN haar 
such that a, Fine) = Pre). 


Equation (93) is differentiated with respect to A and the result is 
ota 2 ' =e 
seh + Bora] > 


where Pp = P (x, 6,7, y) | 


If the limit as H tends to zero is taken, the equation reduces to 


IE a =/2 a 
abr = Q Ara) zer, Pexe - a Bix 18 ia 


The solution for FRE T) is chosen to be 


E, (%, ET) = UlxE) exp [a Eye] œ% 


Hence the stability of the system may be analyzed by inquiring whether 
A (€) is negative or positive. The value of ACE) is described 


by a Taylor series expansion about the parameter E and is given by 


AEE ER (97) 
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If the proposed solution for T (HT, €) (.e., Eq. (96) is sūbsti- 


wed into Eq. (95), then 


an A 
TUCK, €) T EZ mek (Sees ey Pix e) (98) 
3 | 
a Lex | Ulxe)=O ind 
U (xé) = 0O - on SD 


Taking the limit as E approaches zero results in the expression: 


2 O a ia 
V U (xo) (Rs de ze 2.65 Uio) on 
Uu Cx o) = (© on dD. 


A solution to this equation is offered by a modification to the R. Courant 
and D. Hilbert perturbation method (provided in Appendix B). Their 


approximated solution is given by the series 


BR CU Git — 100 
ee CC) tee eee) + (3) Tipe ae 


a e 
Since the contribution of ( 72.) is very small,only the first and 
second terms of the series are retained. The technique provides that 


Ge (X) is the solution to 
T a =- ; 
V cx) F er -x.,) Lux) =O WD 


LA Ox) = © On SD 


> 
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Since the value of A om has already been determined in Eq. (72), 


the only nontrivial solution for U ( XI is 


where Wa is defined by Eq. (73). (Proof of this solution is pro- 
vided in Appendix C). The eigenvalues corresponding to WA mml X) are 
Aoa = A SAN Aa Biene loresche first term tor the series for 


x (€) is defined as 


ETA = De - Aom. (101) 


The linear term of the series of Eq. (100) is given by 





OS 
r=) Kyo 
F oe 
9=135: 


zm 
where Am; =f Be p n ax 


The solution for LA may readily be seen to be 


U = Di | oe +8) 32 Dom = a Vrn +] 


eee 


te ay aan Ym ce 


This solution may be put into the form afforded by Eq. (81) 





to = ee A (102) 
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If Eq. (101) is considered, some conclusions as to the stability 
of the reactor for very small € canbe made. For the fundamental 


mode of the reactor (M= 4 ) the value of Kom is 


K oim = AR, = Wise 
(I) U-m). 


The possible values of Bos oe ake 


i 


AS = KO pn =1 (102a) 


= -27° m = 3,57... 


Thus fer values of E near zero the first term of the series for 
ol (€) indicates that the flux is neutrally stable for the fundamental 
mode of the deviation A 1 ) and stable for the higher modes of the 
deviation (W741). The question of neutral stability obliges an in- 
vestigation into the remaining terms of A CE). 

some conclusive results may be obtained concerning the higher 
modes of reactor operation. The possible values of AE) for the case 


where M > 1 are 


2 
A (tE) = “er pm > mM 
= a | fran = N (102b) 
Z 
= 29 ar 


02 


Since ol (E) may be negative, the higher modes are unstable. 
The linear term in the series expansion of (E) is evaluated 


by differentiating Eq. (98) with respect to € and thus 


; oa Yi 
Vx €) +16) = zaPfıye OT) 


- de) uno + [6 - Epi) - HE] 


*U(%,6) = O. 
The value of A(O) has already been determined to be zero. If the 


limit as E approaches zero is taken, then 


ae A mi S 
Y LOCO) SE e PAZ P, = Kom + Lz €X,0) (103) 
Quake F 
=(¢ P exo) Pixa + &) ULXO), 


The solvability condition requires that 
=] 
Q. /2 . 4 Ty (104) 
4 P ox) Pox) + &) UA Ain) aX =O 
D 
e H 
where LA is the homogeneous solution to Eq. (103). If the value 
of a ìs substituted into the homogeneous equation, the result 
is an equation similar to Eq. (99), thus 
MH eh 
Um 7 U m I 


- K 
When U en and WA m are used to evaluate the solvability 


$ 


condition the value of A is developed. The details of this procedure 


are developed by using the results for P, cx) ; p ( x) ang 


73 





{i (x) obtained from Eqs. (74), (85), and (102), respectively 


. H 
and lA (x)  . Substituting these functions into Eq. (104) gives 


- Ya Y, . _ 
J G An BN + o) EL UT. dx =o 

i -Y, Y 2 
af 7. dx A SE dx = O. 005 


Consider the coefficient of CX in view of the definition of JA 


O 


given by Eq. (81) 








= 1T OCR 
— 


An approximation for the linear term of Eq. (105) may similarly be 


obtained to be: 
E, 2 LO Y 
iel We je ly ax o Be fy, io dx 
D 
+0 (01) 


if terms of the second power of ( a ) and larger are considered negligible. 
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iaiese approximations are substituted into Bq. (105) and 
° “fa, 
gs -taly Y? dy (107) 
i D MA > 


The stability of the system is determined from the evaluation of 
by substituting the results of Eqs. (101) and (107) into the expansion 


OLEG. (97). Hence, 


2 Q % 
6 =(%) (m’- m’) - a RJ a Pn da tee 


e 


The additional term enables the further discussion of the stability 


of the fundamental mode of reactor operation. Thus, for M = í. 
ic we 2 
_(¥, a ty 1, A | a 
ato =(Y%) (tm) - €) NV dr te", 


For all values of AMM >Å the first term becomes dominant and forces 
the value of Œ(€) to be negative. However, the choice of NAL 


defines the value of X (€) ; 
$ 
OS y Re 
K€) = O-4e MI To aX « se 
This integral may be evaluated and the result is, 
dé) = -€ pr (279) (111) 
JA, (279) . 


This result provides the following conclusions; 
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a) negative reactivity 


oA (€) < O for positive € 


o (€) > O for negative € 


These results indicate that the equilibrium solution described by 
an increase in the flux over the initial flux distribution will be stable. 
On the other hand, the other equilibrium solution, which was determined 
to be unacceptable, is unstable. 


b) positive reactivity 


04 (€) > Ô for positive € 


o (E) — O for negative € 


The equilibrium state described by a negative value of € 
is observed to be stable. The other equilibrium state which was con- 
cluded to be unrealistic is unstable. 

In summary, the stability of a nonlinear reactor was studied by 
solving the boundary value problem and then considering the effect of 
an initial deviation on these solutions, called the equilibrium states. 
The possible equilibrium states are those described by Eq. (91) with a 


unique eigenvalue for each state given by Eq. (92). However, the 
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stability criterion predicts that only one of these states are stable and 
that state will be physically observed in the reactor. For a negative 
temperature coefficient the stable state is one described by a positive 
€ . For the case of positive temperature coefficient the negative 
value of € represents the stable condition. Other equilibrium states 


are unrealistic. 


©. EQUILIBRIUM STATES FOR A FINITE CYLINDRICAL REACTOR 

The analysis of the equilibrium states of the finite cylindrical 
reactor can be done in a manner similar to that done for the slab reactor. 
The equilibrium states and the associated eigenvalues will correspond 
to those developed for the slab reactor. However differences will be 
noted because of the differences in the geometry. The flux distribution 
depends upon a radial coordinate, C , which varies from the center 
line to the outer radius, (3 , and a longitudinal coordinate, Zz : 
“woach varies from -l to +l. 

The equilibrium solution is developed from the governing equations 


CRT) 
which are simplified by requiring o Pe ad Von = © A 


The governing equations, Eqs. (3) and (8), reduce to the form, (112) 


2 = = 
V Perz) +[\+ O( Per ES) Pie) = O 
in D 
Pirz) =O on Ay. 
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The perturbation solution is evaluated by a Taylor series expansion : 
about an initial reactor state, Pp (v,2) and NE . If € is the 
perturbation parameter, which will subsequently be defined, the resultant 
equilibrium state is 

o i 2 . 0 
Pre) = P (G2) + € Paz + 2€ Per a) +... (13) 
e 4 2 se 
Ce) mer EN ED ENTE 


de 


The initial states, Der 2) , of the reqctor is defined as the solution 


to the system: 


€ 
FT "Luz + A Pr =0 ape) 
Pina) = O Scorn 
The function Pira) is determined in a manner similar to that already 


meewviced for Eq. (43). The result is 


Pies = An E As UA 


If the shape function YW, is defined for the cylinder as 


Vy, = Jolen) Cos ("Ja 2) al 
then rz) = An Y : (116) 


The eigenvalues are also determined by 
E ee es 
MAT 
1: (Ava) +( In) (117) 


th 
where Pm is the n root of the zeroth order Bessel function. 
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The linear terms of the expansion of Eq. (113) are determined by 
differentiating Eq. (112) with respect to € -© and then taking the limit 


as € tends to zero. The result of these operations is 


. OL =V. > 
TP (ye) + (>. - F¥ wa) Pega) = ma 
ay A Pre) in D 


Pina) = 0 on )., 


In order that a solution is obtained for Eq. (118) the solvability con- 


dition (given in Appendix A) requires 


: sH 
J -J P irz) P a) Lyr dr daz2 = 0 (119) 
D 


> H 
where p cr 2) is the homogeneous solution to Eq. (118). The 
>e H 
equation defining P OV,2@) Eq. (120), is solved by a modification 
to the R. Courant and D. Hilbert perturbation method (provided in 


Appendix B). Thus, 
-© H a a 
2 a 2 
G Yina + (A. al Pa) Piz) =O (20) 
; D 
in 


> H 
Pina) = O on SD 


The approximated solution is given by the series 


MN Q | 
(Pro) = U m 01,2) F 2 (Vin (V2) +- Oa . (121) 
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Only the first two terms of the series are retained since the contribution 
2 
of Q is small. The solution for Uns Z) is obtained from the 


equation 


€ . 
NEO Ze le le zer m 


Um (1 7) = O On SD), 


Since the solution for this equation was already determined for Eq. (114), 


the result is similar and 


Uml?) = Eu A 2 


The value of IR is determined by the orthonormal condition im- 


posed on the Un (r, z) UnNcHOn ATUS: 


[u Aye yy oe es = is 
D 


The integral is evaluated and the result is that 
L = A == n , 
x = 2 
¡A M Y ES J, (um) (122) 


The second term of the expansion is defined as 
GO 








mz) = 
Na ) | ee me Ao; Us 
nas 
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where A 
M 


git 
j J. ES e (178) Aj (a) Ey rar az 


ARAYA MT). 


and Vom 


The value of Alm (Y, 2) is 


cS 


rn Cee) E ). ne Ao; Uy Y; 
ghee 








where 
la ~/2 


The solution for p Cr, 2) may be determined by evaluating the 


series of Eq. (121). Thus, 





OS 
= : ay dmj o 
Per z) u Ba O, Tig er O) Y; 
313,50": 
Em 
If the shape function pi is defined as 
M 
= O, pe "i s Ron- Ao ro} OF y; To... (124) 


e 
IPN 


where De, f Nem , Ami , and Nt. are w by 


Egs (122), 17); (123) and (TIS ie aaa O Dex Z) is 
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er H 
PY (r2) = Bm 7: a > 


The values of g CY, 2) mGivenue bq ait), and Eq. (125) 


are used to evaluate the solvability condition, Eq. (119). Thus, 


J-> An% Ba Ta 2dr dz 


i 


O 


or 


A Y, T, r de dz = 0. 


This equation is satisfied by the condition 


NS (126) 


wan this result is substituted into Eq. (118) the solution for Gs t) 
Cne homogeneous solution to Eq. (118) or the solution to Eq. (120). 


Thus the value of Pir Ei 


Pira) z pr “En . (127) 


The value of the perturbation parameter, E , is chosen to represent 
the average magnitude of the linear term in the expansion of Eq. (113) 
weighted by the T function and measured in units of the maximum 


initial flux. Thus 


E = ah ale | Peixe) =P (ox) | ax. (128) 
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If Eq. (128) is differentiated with respect to E and the limit is 


taken as É approaches zero, the result is 


l e ¿yr dr dz 


Am; 2 
-/[e. A, + Y Dan rej O, Y | tyrdradz, 


at 





2 
If the quantity within the bracket is expanded and terms of Ch are 


considered small, then 


2 2 
== ~ 2 dr dz 
[ore a 


= 

T ajon MA E hos O; N, cy r drdz 
41,5," 
in 


O 


Because oí the orthogonality oí the ab, function the second term 


vanishes and 


oe = [0% Y; ¿yr dr dz 
= Í, 


| 


Thus; An = Bm and Eq. (127) may be rewritten as 
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Diz) = Ar: (129) 


The quadratic terms of the expansion for ARS and A (€) 
given by Eq. (113) are evaluated by differentiating Eq. (112) twice with 


respect to E and then let E tend to zero. The result is 


ig y N 
Y Pira) de (>, E > P erz) Pen) = (130) 


y ANS 
- A Ri) Te 2 Ana Pin) in D 


P (1?) = O on SD. 


The solvability condition requires that 
-3 
oe Q. 72 wack . H dor 
í > Pia) + 7 hd) Per gy) Y (2) er draz =0 
D 


PR 
where PL 2) is the homegeneous solution to Eq. (130). This 


was previously determined for Eq. (120) and from that result 


1 Y 


H 
(7,2) = Cm EN. (132) 


lt 


The value of A may be determined by substituting the correct 


functions determined above. Thus, 
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N IR Ain Yon az ig dy r dr dz = 
a f A y ie CAT 2yr dr dz. 


The coefficient of A is evaluated as follows: 


Er rn, 271 dr dz 


2 A mf hlo it 2) 3 Y | 
3156, 
x ¿nr ar Az 
z Am 
— 


where LAN is defined by Eq. (122). The integral on the right hand 


side is evaluated as follows: 


) AS 
TARA Ta 2yr dr dz 


mE a „NY 3/2 3 
= ca O, 2yr dr dz 


a e 6 
if terms multiplied by Q are considered small. The value of A 


= 3 y 
\ = m] veo, ¿rr dr dz. (133) 


is 
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For the fundamental mode (M=1 ) 
ee 56 3 A 
= ZA 37 € (134) 
A Yo VA 


The various terms of the expansion of Eq. (113) may be combined 
to give the equilibrium states of the reactor. Thus, substituting Eqs. 


(116) and (129) into the expansion for P (r, zZ, €) gives 


P(r, 2 €) = Am Ww + E fhe + Olé | dos 


and substituting Eqs. (117), (126), and (133) into the expansion for 


A(E) gives 
2 | 
re) = (MA) + (TR) 
I 
rte A 2 Jn On 2yr dr dè +e. 


For the fundamental mode of the reactor (M=J_ ) the integral is 


evaluated and 
2 2 2 2.563 ar (136) 
Ne) (Fr) a) HE “ee fart 


It is apparent that for any mode of reactor operation there will be two 
equilibrium states defined by Eq. (135) and each of these states will have 
unique eigenvalues defined by Eq. (136). However, only one of these 


values will represent an acceptable reactor state. 
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1D., STABILIBTANALYSSSOFRZTHE FINITE CYENDRICALCREACTOR 

The stability criterion for a finite cylindrical reactor with negligible 
fraction of delayed neutrons is determined by the methods previously 
considered. An initial deviation, YN , due to temperature feedback 
is introduced into the equilibrium states of the reactor and the temporal 
effect is noted. 

The governing equation has the form given by Eq. (93) where the 
functional dependence of the variable is changed to accommodate the 


cylindrical geometry 


one 7) 


(137) 


a 
E =- J KE A - al Pis) - Le) Rca 
ino D 
Pinar) = O on S&P 


P (1,2,0) = A (1?) . 


The deviation is described by a Taylor series expansion about the 


equilibrium states, Eq. (135), given by 
(138) 


2 
Paren) = Penze +y leze) +z% Rauen) +0 


such that 


= P (reeg) = Py2,¢) 
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Eq. (137) is differentiated with respect to u and then the limit is 
taken as the value of H tends to zero. The result of these operations 


yields the equation 


2 
O fy (RET) = Y E, (G3,€,7) +|) (139) 
e 


Ye 
= = Pire er) or. (ea 442, 6,0) In D 
, (12,67) = O De 


A form for the solution of A (CLET) is Chosen to be 
h Fuer) = uU Ghe) exp | xe) a : a 


The stability of the system may be determined by investigating the sign 
of xX (E) . If the sign of oA (E) were positive, then the exponential 
would increase with time. Hence, the contribution of Pa in the 
series of Eq. (138) would be unbounded with time and the system would 
be unstable. On the other hand, negative value of xX (E) wouid indicate 
that the Py contribution would die and stability would be assured. 


The value of A (€) will be determined by a series expansion given by 
> | Le 
= - pena ee 9 (141) 
KE) = Ll, t+ OX +7 A+ 


The following will be directed toward evaluating each of the terms of 


this expansion. 
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If the proposed solution for A lr eee €) pec el Oe as 


substituted intoEq. (139) the result is, 


E a 
Vu t E =~- Kle) +t F rre) (142) 
el D |u = O wD 
Beeren = O on SD 


where LL = UL Cv E €) . The first term of the series of A (€) 


is evaluated by taking the limit of Eq. (142) as © approaches zero. 


Thus, 
i —/2 
Vu t EAN = Os > tat a] (143) 
LTE) = O in D 


The solution to this equation may be determined by the approximation 
given by the R. Courant and D. Hilbert perturbation technique. A 


solution to Eq. (143) is developed about the series 


aa OL 2 
Una) = UG + Z Mm + 0O(0). 


The solution for UT 2) is determined about the solution to the 


equation 


Vo hia) + Era - 0.0) Giray = O ind 


LE) 


i 
O 
O 
> 
on 
U 
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It may be noted that the solution to this equation is similar to that ob- 

served in Eq. (114). However, by the results of Appendix C one may 
conclude that the solution for U (r 2) is given by Ur 2) = Ur 
where the function A is defined in Eq. (115) and Ne com 


Thus, 


Ko = Aoem 7 Aom (144) 


and the linear term of the NE) series of Eq. (141) is defined. (A 
discussion of Xeo will subsequently be made.) Since the La) 


function must be orthonormal the constant U zn may be evaluated 


by the normality condition 


f Tn Ya nr ar de =1. 
D 


a 
The value of U m is 


7 * = UL Bei (122) 
Ulim = um A A 


a | 
~ 
-5 
eN 


which has previously been determined by Eq. (122). Thus 


is 


Inne) = On > 


The second term of the expansion \AW,2) is defined as 
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ai Y eae 
313,5": J 


# me 
or 


co 
Nor we =) Nom — roj 9, Y, 
(¿NES "> 
ie | 
where dm; is defined by Eq. (123). The solution for the U (62) 


function from Eq. (143) is determined as 


= E 
uta = E a a2 Na 





I 
G 
3 
=== 
= 
i 
+ 
NJM 
me 
xio 
3 3 
E 
ce 
Lae 


= Um Im (145) 


where T isade ined by La. (124). 
The value of of omrmm as evaluated from Eq. (144) is rewritten 


below 
EN = eke = Nom (146) 


If the value of NE as given by Eq. (117) are used, the result is 
| / 2 2 | y 2 a 
_ [— A e = (1463) 
OC = (7) Hm Hem u Y, (m ORJ 
For the fundamental mode ( M = 1) this equation reduces to the form 
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Kam = (rn) (24048 - Um) ATA) e il). 


For values of É near zero the values of PAA are 
Xeon = O m m 


<-2 mM >M : 


Thus for all modes of the deviation except the fundamental the sign of 
NS Fin is negative and any initial deviation of the system will 
decay with time. However, for the value MM=M the first term evalua- 
tion of the stability indicates that there is a neutral stability for the 
fundamental mode. Thus it is necessary to investigate the higher order 
terms of the expansion of A (€) . For the higher modes of reactor 


Operation oO (e) may have the following values, 


2 
A (ed Ze) m im 
= mM EM 
Z 
< -lT IM > We. 


Thus the value of olle) may be positive when MN for M”>Í, 
Since the value of Coan m is so large when compared to the linear 
term (order of E ) of ol (€) , one may conclude that the first 
term of the series is dominant and the higher modes are unstable. 
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The linear term in the expansion of Ale) Eee) ana 
be evaluated by differentiating Eq. (142) with respect to E and taking 


the limit as € temas to zero. Thus, 


=e ia (147) 
Ä F ! 
(Y O Poner Perz) + o U) 
mo D 
Lk (42) =D on 6D | 


The solvability condition requires that 


z 
> dl ~% = f H 
Fé + ¿0 Pa) Pira) tee) une) Zyrdrd® =O 
> 


Y 
8 
where (4 (1,2) ¡Stereo nesencous solution to Eq. (147). This 


solution was previously developed for Eq. (143) and the solution is 


LL (42) = Den Ea 


Y 
With the solution for UL (12) and the folutions for (4. (1,7) 
P, (1,2), and Por, =) gien y Eqs. (125), (116), and (129) 


the value of X may be evaluated from the solvability condition. 


Thus, 
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. i 2 | 
Um ml Te 2yrdr dz = 
A : Y 2 
-¢ A Um Emp Y Po Er 2 r dr dz 
a S y 


The coefficient of OX was evaluated for Eq. (128). The result in- 


2 
dicated that if the terms of order Q are considered very small, 


the integral was equal to one. The value of X is 


ua -% 2 
A = 4 A, Y 7. To eyr dr dz, (148) 


Since the fundamental mode of the reactor is of greatest interest, the 


integral of Eq. (148) is evaluated. for MAIZ M = 4 “Thus, 


y = Q ae 3 
Xx ae] Y T 2yr dr da 


Ye 
a | a4, 2yr dr az 


2 373 A 


i 


ile 


are 


The results of Eqs. (146) and (148) may be combined into the expansion 


of OX CE) , Eq. (4D, 
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ELEND (149) 


-3 
Ma fa 2 
et al % Pa Lo, 27 r dr da 


and for the fundamental mode 





END = (4) (24078 - u2) (IN t-mè) 
S75 


a 2 
er OK Nie (150) 





From the previous analysis values of AN > j have been found to 
favor stability. For the value of /M=Í the following becomes 


evident, 





Q- 2 
A(E) = =€ “E Ya + OC), (151) 


For the case of negative reactivity the sign of ol (E) will be deter- 
mined by the value of E . A positive value of a will permit a 
negative sign of XE) . The deviation associated with of (€) 
will decrease with time and stability results. For positive reactivity 
the sign of ( OL ) is changed and the negative value of a permits 
stability. The other states of reactor equilibrium are unrealistic as 


determined by previous arguments. 
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APPENDIX A 


THE SOLVABILITY CONDITION FOR NONHOMOGENEOUS 
BOUNDARY VALUE PROBLEMS 


Given the equation 


TDT Dex) T A Q(x) = Ls in D (A~1) 


? 


Pix) = O on SD 


prove that this equation has a solution if and only if the function fex) 


satisfies the condition: 


AH 
[Fr Yixy AV =O 
D 


where 


2 H H 
V Pox) + APixn =O n D (a-2) 


H 
Pix) = O on SD 


H 
Multiplying Eq. (A-1) by Pix) and integrating throughout the 


domain the following results: 


H 2 H 
| (ten V Pix) f » Poy Lex) ) AV (A-3) 
D 
- | Dex fon ay = 
D 
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2. H H 
Brom Eq. (A=2) JZ YP (x) = -A LPi X) and this is substituted 


into Eq. (A-3). Thus, 


H 2 
| (ec V Pix) - Pix) Y Pex) dV i 
D 
-| Dix) Pex) av. 


By Green's Theorem the left hand side of Eq. (A-4) may be transformed 


into a surface integral: 


E H 
[nom ven em vet) as 


=f fox Pox dv. 


Thus, if Pix) is to be the solution to Eq. (A-1), the homogeneous 


boundary condition on SD must be satisfied. Therefore, 


A 
|, Lo) Pix) NA O (A-5) 


O E.D. 
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APPENDIX B 


A PERTURBATION METHOD FOR 
LINEAR EIGENVALUE PROBLEMS [7] 


The solutions to equations of the form 
(BS 
2 = Pere aa . 


Um = O on & D 


where E v(x) represents a small perturbation , ( 





Er co] E Am), 


are provided by a series expansion about the solution of the equation 
È : 
Ne ny ene ec: O in D (B-2) 


Um = O on &D 


The eigenvalues are held constant while the eigenfunctions are perturbed 


in accordance with the Taylor series given by 


- mM do 
Um = Um + ENm + TE Um a 


Hane series is substituted into ra ap 1) the result is 
2 \ 2 
Ve ee E I.) 
2. 


- Erım (um + Ena + ¿CE Ut: 


+ rn (Um + E Nm +3 € mt) 
ZAG. 


98 


By separating terms of the same power in E , the following equations 


are obtained: 


VA O (B-4) 

et: 

Y Na ir N Ma NEL) U (B-5) 
Z 

V Un t Am Um = VY(x) Um.  (B-6) 


Multiply Eq. (B-5) by LA 2 and integrate over the volume, av 


flus Vz t Am Lp Nin) av of roo uu AV 
D D 


fu, T - Nm Ts Zee V us) AV (B-7) 
> 


Anm AN =f V(X) Un Ue dV. 
D 


By Green's Theorem the following volume integral may be transformed 


into a surface integral 


| (Team = Nm vu,) AV 


ii 


[0 (Ua Trim m vu dS 
5 


However, the boundary conditions specify that Up and Nr 


are everywhere zero on the surface and the surface integral vanishes. 
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Equation (B-7) may be rewritten as 


[er YU + AmUz Ma) dV = dm 
where 

O = J DEE U Y. (B-8) 
Equation (B-2) reveals 


2 
Y uy = ee, 


which is substituted into Eq. (B-8) to give 


de = RE Om = Ch ms (B-9) 


where 
Om: = A Us dV 5 
D 
The value of the integral Q MAN is evaluated if the 
L 
normalizing condition J, A AV = lk is imposed. It is 
clear that 


(Un + E Na + e oo) AV =]. (BELO) 
D 


By equating powers of É on the right and left sides of Eq. (B-10) 


f únax =1 


[ um m dx = © 





- 


Or ewer =0. 


The function Mm (X) is determined by a series of orthogonal 


functions given by 


od 
CR) e Da; U, y (B-11) 


Multiplying both sides by us and integrating over the domain 


[an u, qY =| È bmj uj Uo ave 
D j= 


Assuming that the expansion exists and that the integral of the infinite 


series is equal to the sum of the integrals, 
Od 
Ch mg =, De Í u dV 
iat 3 D 


Amm =O, 


Since 


f u; Umay = © Ge 
li; Up AV =i =X 
and f. J x J ? 
the constants Dmi are evaluated as Dmi = Q mj 3 
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nz Given 
D AmA) 


Thus the series of orthogonal function given by Eq. (B-11) may be 


However, from Eq. (B-9) it is obvious that b 


evaluated as 22 


Ginna 
A A ( 
Mm OR Am - Aj u, x) 
y=1?,7-- 
Am 


The series of Eq. (B-3) may be evaluated as 


er < mj | 
Br e) Am Ay Uy ex) +. 


ia 
where Am; = [roou. U4 dy . 
D 
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APPENDIX C 


Given the equation 


LV T A Vn (x) = © m D (C-1) 


Vin A) 10 on sD 


with eigenfunctions D(x) and eigenvalues we - Prove 


that if 


2 
Y Pix) A Pox) = O m D (C-2) 


aa on Da 


where A Æ een is a new physical parameter, then the nontrivial 


solutions to Eq. (C-2) are 


Pix) = Pry OX) Lie Ba 


where fm and /M are independent subscripts. 


Multiply Eq. (C-2) by W (xX) and integrate throughout 


the domain. 


/ 7 (x) T Pix) dV + af Loo Pix) AY =O 


z 
By adding and subtracting Pio Z Y, exe) 


KA V Poo - Poo Van) dv e 


J, 2e Z Vow av +f Yo Ax avy =O. 
D 


From Green's Theorem the first integral may be transformed into a surface 


integral by 


Jo T Pix -= Pin T) dv 


= fn (Yr (x) J Px) — Wx) I Yn (x) ar. 
S 


However, from the homogeneous boundary conditions of A, (x) = Lox) =O 


on SD the surface integral vanishes. Thus Eq. (C-3) is 


[e A (x) dV + ben Dix) ay so 


2 
From Eq. (C-1) V Wy, (x) in DR Vn Cx) j 
therefore 


u An | Poo hos ay + rf W(x) Lx) dV = 0 
D 
D 


If (A 5 a) F O , then the integrai must vanish throughout 
the domain and W(x) is orthogonal to W,, (X) . Hence, 
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Dix) — vy, (X) 


and the eigenvalues of Eq. (C-2) (i.e., A ) are 


NI ns 
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